Analysis of Discrete and Hybrid Stochastic Systems 
by Nonlinear Contraction Theory 
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Abstract — We investigate tlie stability properties of discrete 
and hybrid stochastic nonlinear dynamical systems. More pre- 
cisely, we extend the stochastic contraction theorems (which were 
formulated for continuous systems) to the case of discrete and 
hybrid resetting systems. In particular, we show that the mean 
square distance between any two trajectories of a discrete (or 
hybrid resetting) contracting stochastic system is upper-bounded 
by a constant after exponential transients. Using these results, we 
study the synchronization of noisy nonlinear oscillators coupled 
by discrete noisy interactions. 

I. Introduction 

Contraction theory is a set of relatively recent tools that 
provide a systematic approach to the stability analysis of a 
large class of nonlinear dynamical systems [1], [2], [3], [4]. A 
nonlinear nonautonomous system x = f (x, t) is contracting if 
the symmetric part of the Jacobian matrix of f is uniformly 
negative definite in some metric. Using elementary fluid dy- 
namics techniques, it can be shown that contracting systems 
are incrementally stable, that is, any two system trajectories 
exponentially converge to each other [1]. 

From a practical viewpoint, contraction theory has been 
successfully applied to a number of important problems, such 
as mechanical observers and controllers design [5], chemical 
processes control [6], synchronization analysis [2], [7] or 
biological systems modelling [8]. 

Recently, contraction analysis has been extended to the case 
of stochastic dynamical systems governed by Ito differential 
equations [4] . In parallel, hybrid versions of contraction theory 
have also been developped [3]. A hybrid system is charac- 
terized by a continuous evolution of the system's state, and 
intermittent discrete transitions. Such systems are pervasive 
in both artificial (e.g. analog physical processes controlled 
by digital devices) and natural (e.g. spiking neurons with 
subthreshold dynamics) environments. 

This paper benefits from these recent developments, and 
provides an exponential stability result for discrete and hybrid 
systems governed by stochastic difference and differential 
equations. More precisely, we prove in section and |III] that 
the mean square distance between any two trajectories of a 
discrete (respectively hybrid resetting) stochastic contracting 
system is upper-bounded by a constant after exponential tran- 
sients. This bound can be expressed as function of the noise 
intensities and the contraction rates of the noise-free systems. 
In section |IV] we briefly discuss a number of theoretical 



issues regarding our analysis. In section [V] we study, using 
the previously developped tools, the synchronization of noisy 
nonlinear oscillators that interact by discrete noisy couplings. 
Finally, some future directions of research are indicated in 
section |VT] 

Notations The symmetric part of a matrix A is defined 
as As — \ (a + A'^). For a symmetric matrix A, A,nin(A) 
and A,nax(A) denote respectively the smallest and the largest 
eigenvalue of A. A set of symmetric matrices (Ai)^^^ is 
uniformly positive definite if 3a > 0, Vi G /, Amin(Ai) > a. 
Finally, for a process x(t), we note Ex(-) = E(-|x(0) = x). 



II. Discrete systems 

We first prove a lemma that makes explicit the initial 
"discrete contraction" proof (see section 5 of [1]). Note that a 
similar proof for continuous systems can be found in [9]. 

Lemma 1 (and definition): Consider two metrics Mj = 
defined over M"' (i = 1,2) and a smooth function 
f : M"2. The generalized Jacobian of f in the metrics 

(Mi,M2) is defined by 

ox 

Assume now that f is contracting in the metrics (Mi, M2) 
with rate /3 (0 < /9 < 1), i.e. 



Vxe 



A„,ax(F(x)^F(x)) < /3 



Then for all u, v G M", one has 

dM.(f(u),f(v))2 </3dMi(u,v)2 

where c^m denotes the distance associated with the metric M 
(the distance between two points is defined by the infimum of 
the lengths in the metric M of all continuously differentiable 
curves connecting these points). 

Proof Consider a C^ curve 7 : [0, 1] — > R"^ that connects 
u and V (i.e. 7(0) = u and 7(1) = v). The Mi-length of such 
a curve is given by 




Since f is a smooth function, £(7) is also a curve, with 
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The chain rule next implies that 
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which leads to 
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= /o V(lSV)F^F(0,|l)du 



(1) 



Choose now a sequence of curves (7„)„gN such that 
lim„^oo -^Mi(7n) = rfMi(u, w). From ([B, one has Vn e 
N, LM2(f(7n)) < vT^-^Mi (7n)- By definition of distance, one 
then has Vn G N, divis (f f (w)) < VPLm^iln)- Finally, 
by letting n go to infinity in the last inequality, one obtains 
the desired result. □ 

Theorem 1 (Discrete stochastic contraction): Consider the 
stochastic difference equation 

afe+i = f (afc, /c) + a{a.k,k)wk+i (2) 

where f is a M" x N ^ M" function, cr is a M" x N ^ M"'' 
matrix-valued function and {wk, fc = 1, 2, . . .} is a sequence of 
independent d-dimensional Gaussian noise vectors, with Wk ^ 
-^(0,Qfc). 

Assume that the system verifies the following two hypothe- 
ses 

(HI) the dynamics f (a, fc) is contracting in the metrics 
(Mfe, Mi,+i), with contraction rate /3 (0 < /3 < 
1), and the metrics (M/c)fegN are uniformly positive 
definite. 

(H2) the impact of noise is uniformly upper-bounded by a 
constant ^/C in the metrics 



Va, fc (iMfc(f(a, fc), f(a, fc) +(T(a, fc)wfe) < VC 

Let afe and b^. be two trajectories whose initial conditions 
are given by a probability distribution p(xo) — p{ao, bp). Then 
for all fc > 



(dMfc(afc,bfc)) < 
2VC ' 



(iMo(a,b) 



2VC 

+ 

dp{a, b) 



(3) 



where [•]+ ~ max(0, •). 

This implies in particular that for all fc > 

2VC 



E(dMfc(aA;,bfe)) < 



v//3 E(dMo(ao,bo)) (4) 



Proof Let x = (a, b)^ € M^". We have by the triangle in- 
equality (to avoid long formulas, we drop the second argument 
of f and cr in the following calculations) 

rfMfc+i(afc+i,bfc+i) < dMfc+i(f(afc),f(bfe)) 

+ rfMfc+i(f(afc),f(aA;) + cr(afc)wfe+i) 
+ (iMfc+i(f(bfe),f(bfc) + cr(bfe)w;;^_^i) 

Let us examine the conditional expectations of the three 
terms of the right hand side 

• From (HI) and lemma [T| one has 

IEx(rfM,+i(f(afc),f(bfc))) < v/^Ex(dM,(afe,bfe)) 
. Next, from (H2) 



IEx(rfMfc+i(f(afe),f(afe) + (T{a.k)wk+i)) <VC 

and similarly for dM^+i (f (b^), f (b^) + a{hk)w'^._^_^). 

If we now set u^ = Ex((iMfc (afc, b^)) then the above 
implies 

Uk+i < \/Puk + 2VC (5) 

Define next = Uk~ 2\fC l(\ — \fP). Then replacing Uk 
by Vk + 2\/C/(l - Vi?) in © yields 

This imphes that Vfc > 0, Vk < voy/p'' < [vo]+y/p''. 
Replacing Vk by its expression in terms of Uk then yields 

+ 



Vfc > life < 



2VC 



uo ■ 



2VC 

1- V? 



which is the desired result. 

Next, integrating the last inequality with respect to x leads 
to (O. Finally, (|4|l follows from (|3]l by remarking that 



dMo(a, b) - 



C 



dp{a, b) < 



y dMo(a, b)dp(a,b) = E(dMo(a,b)) □ 

Remark In the particular context of state-independent met- 
rics, hypothesis (H2) is equivalent to the following simpler 
condition 

Va, fc tr (cr(a, fc)'^Mfe+icr(a, fc)Qfc) < C 
Also, for state-independent metrics, one has 
rfMfc(afc,bfc)2 = ||afc - bfell^^, = (a/c - bfe)^Mfe(afc - b^) 
which leads to the following stronger result instead of (|4|i 

2C 



El 



|afc-b,||2,J <-_-^+/3^-E(||ao-b, 
III. Hybrid systems 



'OllMo) 



We have derived above the discrete stochastic contraction 
theorem for time- and state-dependent metrics, contrary to 
the context of continuous systems, where the state-dependent- 
metrics version of the contraction theorem is still unproved [4] . 



We now address the case of hybrid systems, but due to the 
current Hmitations of continuous stochastic contraction, only 
state-independent metrics will be considered. 

For clarity, we assume in this paper constant dwell-times, 
although more elaborate conditions regarding dwell-times can 
be adapted from [3]. 

Consider the hybrid resetting stochastic dynamical system 

V/c > a(fcr+) = fd(a(fcT-), k) + crd(a(fcr~), k)wk (6) 

yk > 0, yt e]kT, (fc + l)r[ da ^ fc(a, t)dt + ada, t)dW 

(7) 

All the contraction properties below will be stated with 
respect to a uniformly positive definite time-varying metric 
M(t) = @{t)^ @{t). Furthermore, it will be assumed that for 
all fc > 0, M is continuously differentiable in Jfcr, [k + 1)t[. 
Finally, M(fcr^) and M(/ct+) will respectively denote the left 
and right limits of M(i) at t = kr (and similarly for 0). 



A. The discrete and continuous parts are both contracting 

Theorem 2 (Hybrid stochastic contraction): Assume the 
following conditions 

(i) For all fc, the discrete part is stochastically contract- 
ing at kr with rate (3 <1 and bound Cd, i-e. 

Va e M" A„,ax (F(fc"r)^F(fc"r)) < (3 

where F(fcT) = &{kT+)^{a,k)®{kT-), and 

Va e M" tr ((7d(a, fc)^M(fcr+)(7d(a, fc)Qfc) < Cd 

(ii) For all fc, the continuous part is stochastically con- 
tracting in ]fcr, (fc + 1)t[ with rate A > and bound 

Cc, i.e. Va G R", Vt e]fcT, (fc + 1)t[, 



A, 



|0(t) + 0(o|)0-(O) <-A 



ti{a,{a,tfM{t)<j,{a,t)) < Cc 



(8) 



Let a{t) and h{t) be two trajectories whose initial con- 
ditions are given by a probability distribution p(x(0)) = 
p(a(0),b(0)). Then for alU > 



En|a(t)-b(t)||^(, 



< 



Ci+E(|!a(0)-b(0)||2,(o))/3L*/r 

where C, = ^^^"tillgil^f.y^^^- and n = pe'^'^ . 

Proof For all t > 0, let u{t) = E (^\\a{t) - h{t)\\l^^^^^ and 
let us study the evolution of u{t) between fcr+ and (fc+ 1)t+. 

Condition (ii) and theorem 2 of [4] yield 

ii((fc + l)r-) < ^+u(fcr+)e-2A^ (9) 
Next, condition (i) and theorem [1] above yield 



Substituting ^ into leads to 



u((fc + l)r^ 



< 



2Cd 

-13 



'it 



■ l3u{kT^ 



-2Xt\ 



- |^ + ^+/3e-2^-u(fcT+) 



Define Z^i = |^ + ^ and = it(fcr+) - Di/(1 - n). 
Then, similarly to the proof of theorem [T] we have Vk+i < 
riVk, and then Vk < [I'o]^. which implies 



u{kT^ 



< 



< 



Di 

i-n 

Di 



Now, for any i > 0, choose fc = [</tJ. Then 

u{t) < ^+^,(fcr+)e-2Mt-M 
A 

A 1 — ?'i 



< 



A 



u(0+)/3'^e 



((fc + l)T+)< f^+/3?i((fc+l)T-) 



(10) 



l-n 

which leads to the desired result after some algebraic manip- 
ulations. □ 

B. Only the discrete part is contracting 

Let us examine now the more interesting case when the 
continuous part is not contracting, more precisely when A < 
in (l8|. For this, we shall need to revisit the proof of theorem 2 
in [4]. 

Theorem 3 (Case A = Q): Assume all the hypotheses of 
theorem |2] except that A = in dS]). Then for all i > 

E(||a(f)-b(t)||2,(^)) < 
C2+E(||a(0)-b(0)||^(„))/3L*/-J 

where C, = ^^^+2/3(1-^)^^^ 

Proof As in the proof of theorem 2 in [4], let 

V{x,t) ^V{{a,hf,t) = (a-b)^M(t)(a-b) 

Lemma 1 of [4] is unchanged, yielding (see [4] for more 
details) 

yt e]fcT, (fc + l)r[ Iy(x(t), t) < 2Cc 

where A is the infinitesimal operator associated with the 
process x(t) (see section 2.L2 of [4] or p. 15 of [10] for 
more details). 

By Dynkin's formula [10], one then obtains for all x e M^" 

E,l/(x(<),i)-y(x,fcT+) = E,/* Iy(x(s),s)ds 

< E^jl^2Ccds 
= 2Cc{t - kr) 

Integrating the above inequality with respect to x then yields 

yt e]fcT, (fc + l)r[ u{t) < 2Cc(t - fcr) + u(fcr+) 
In particular, (|9]) becomes 

w((fc + 1)t-) < 2CcT + u{kT+) 



which leads to, after substition into (fTOl l. 

u{{k + l)r+) < + 2/3CcT + l3u[kT+) 



This finally implies 

u{kT+) < 



1-/3 



u(0+)/3'= 



The remainder of the proof can be adapted from that of 
theorem 12] □ 

Theorem 4 (Case A < Oj; Assume all the hypotheses of 
theorem ID except that A < in (|8]l. Let k = [t/T\. There 
are two cases: 

. If /3 < e"2|A|r^ jjjgjj ^2 ^ /3e2|-^l^ < 1. For alU > 



E(^||a(t)-b(t)||^(,)j < 
C3+E(||a(0)-b(0)||^(„))e2|^IV,^ 

Where G3 |A|(i-;3)(i-r2) ■ 

If/3 > e-2|^l^, then there is - in general - no finite bound 



on E 



(l|a(i)-b(t)|| 



M(t) 



as < ^ +00. 



Proof One has now for all t gJ/ct, (fc + 1)t[, 
AV{xit),t) < 2|A|F(x(t),t) + 2Cc 
with |A| > 0. By Dynkin's formula, one has, for all x £ R^" 



ExV"(x(i),t)-\/(x,fcr+) < Ex / i2\X\V {x{s),s) + 2Cc)ds 

Jkr 

Let now g{t) — ExV(x(<), t). The above equation then yields 



g{t) = y(x, /cT+) + 2Cc(t - kT) + 2|A| / 5(s)ds 

Jkr 

Applying the classical Gronwall's lemma [11] to g(t) leads 



to 



g{t) < V"(x,A:T+) + 2Cc(t-fcr)+ 

2|A| Jl {V{x, kT+) + 2C,s) exp (/j 2|A|d^i) 
= W\ (e'l^l^'"''"^ - 1) + ■l^(x, fcr+)e2|^l(*-'=-) 

Integrating the above inequality with respect to x then yields 

Vt e]A:T, (fc + 1)t[, 

|A| V / 

which implies 

u{{k + l)T+) <D2+f3e^^^^^u{kT+) (11) 

where i?2-f^ + ^(e2|^l--l). 
There are three cases: 

. If /3 < e"^!^!"", then r2 = /Se^l^l^ < 1. By the same 
reasoning as in theorem [T] one obtains 



D 



u(kT+) < — + u(0+)r. 

1 - r2 

The remainder of the proof can be adapted from that of 
theorem |2] 



If ^ = e-^l^l^, then (dB reads 

U((fc+1)T+) <D2+u{kT+) 

which implies Vfc > 0, u{kT^) < kD2 + ^(0"*"). From 
this, it is clear that there is - in general - no finite bound 
for u(fcr+). 

If /3 > e^^l^l'^, then r2 = /3e^l-^l'^ > 1. By the same 
reasoning as in theorem [T] one obtains 



u{ki 



D2 



r2 



1 



D2 



7*2 



1 



m 



Since r2 > 1 in this case, it is clear that there is 
general - no finite bound for u(fcr+). □ 

Remarks Theorems [3] and |4] show that it is possible to 
stabilize an unstable system by discrete resettings. If the 
continuous system is indifferent (A = 0), then any sequence 
of uniformly contracting resettings is stabilizing. However, it 
should be noted that the asymptotic bound C2 ^ 00 when 
/3 — > 1. In contrast, if the continuous system is strictly unstable 
(A < 0), then specific contraction rates (depending on the 
dwell-time and the "expansion" rate of the continuous system) 
of the resettings are required. Finally, note that in both cases, 
the asymptotic bounds C2 and C3 are increasing functions of 
the dwell-time r. 

IV. Comments 

A. Modelling issue: distinct driving noise 

In the same spirit as [4], and contrary to previous works on 
the stability of stochastic systems [12], the a and b systems 
considered in sections HI] and |lll] are driven by distinct and 
independent noise processes. This approach enables us to 
study the stability of the system with respect to variations 
in initial conditions and to random perturbations: indeed, two 
trajectories of any real-life system are typically affected by 
distinct realizations of the noise. In addition, this approach 
leads very naturally to nice results regarding the comparison 
of noisy and noise-free trajectories (see section IIV-BI ). which 
are particularly useful in applications (see e.g. section [Vb. 

However, because of the very fact that the two trajectories 
are driven by distinct noise processes, we cannot expect the 
influence of noise to vanish when the two trajectories get very 
close to each other As a consequence, the asymptotic bounds 
2C/(1 — /3) (for discrete systems) and Ci, C2, C3 (for hybrid 
systems) are strictly positive. These bounds are nevertheless 
optimal, in the sense that they can be attained (adapt the 
Ornstein-Uhlenbeck example in section 2.3.1 of [4]). 

B. Noisy and noise-free trajectories 

Instead of considering two noisy trajectories a and b as in 
theorem [U we assume now that a is noisy, while b is noise- 
free. More precisely, for all A; G N 

afe+i = f(afe,fc) + a{aLk,k)wk+i 

hk+i = f(bfc,fc) 



To show the exponential convergence of a and b to each 
other, one can follow the same reasoning as in the proof of 
theorem [T] with C is replaced by C/2. This leads to the 
following result 

Corollary 1: Assume all the hypothesis of theorem [T] and 
consider a noise-free trajectory hk and a noisy trajectory 
whose initial conditions are given by a probability distribution 
p(ao). Then, for all fc e N 

C 



where 



IE(||afc 



2 
Mfc 



< 



(3^ 

Remarks 



2 

I Mo 



c 



1-/?, 



1-/3 
dp{a) 



(12) 



• The above derivation of corollary 1 is only permitted by 
our choice of considering distinct driving noise processes 
for systems a and b (see section ITV-AI ). 

• Based on theorems |2] [3] and [H similar corollaries can be 
obtained for hybrid systems. 

• These corollaries provide a robustness result for con- 
tracting discrete and hybrid systems, in the sense that 
any contracting system is automatically protected against 
noise, as quantified by (fT2] i. This robustness could be 
related to the exponential nature of contraction stability. 

V. Application: oscillator synchronization by 

DISCRETE couplings 

Using the above developped tools, we study in this section 
the synchronization of nonlinear oscillators in presence of 
random perturbations. The novelty here is that the interactions 
between the oscillators occur at discrete time instants, contrary 
to many previous works devoted to synchronization in the 
state-spac^ [14], [7]. 

Specifically, consider the Central Pattern Generator (CPG) 
delivering 27r/3-phase-locked signals of section 5.3 in [7]. This 
CPG consists of a network of three Andronov-Hopf oscillators 
Xj ~ {xi,yi)'^, i — 1,2,3. We construct below a discrete- 
couplings version of this CPG. 

At instants t ~ kr, k £ N, the three oscillators are coupled 
in the following way (assuming noisy measurements) 

Xj(fcr+) = Xi{kT^) 

+ j(R(ic,+iikT-) + ^Wk) -Mkr-)) 

with X4 = xi and 



R 



Between two interaction instants, the oscillators follow the 
uncoupled, noisy, dynamics 

dxi = i{K,)dt + ^dW 
v2 

'Discrete couplings are more frequent in the literature devoted to phase 
oscillators synchronization, where phase reduction techniques are used [13]. 
However, contrary to our approach, these techniques are only applicable in 
the case of weak coupling strenghs and small noise intensities. 



f(x,;) = f 



3 2 
Hi ■^iUi 

x^ + y^- yf - y^xf 



We apply now the projection technique developped in [7], 
[4] . We recommend the reader to refer to these papers for more 
details about the following calculations. 

Consider first the (linear) subspace M of the global state 
space (the global state is defined by x = (xi, X2, X3)"'") where 
the oscillators are 27r/3-phase-locked 

M = {(R2(x),R(x),x)^ : X e 

Let V and U be two orthonormal projections on A^-'- and 
A4 respectively and consider y — Vx. Since the mapping 
is linear, using Ito differentiation rule yields the following 
dynamics for y 

VfceN f{kT+)^gd{y{kT-j)+j^wk (13) 

v2 

\fte]kT,{k + l)T[ df ^ gc{y)dt + ^dW (14) 

v2 



with 



where 



i(y) = VLx = VL(V'^y + U'^Ux) ^ VLV'^y 
ge(y) = Vf (V^y + U^US) 



(1 - 7)12 7R. 
L = I (1 - 7)12 7R 

7R (1 - 7)12 

f(x) = (f(xi),f(x2),f(x3))^ 

Remark that gd(0) — and gc(0) (the last equality 
holds because of the symmetry of f: Vx, f(Rx) = R(f(x))). 
Thus, is a particular solution to the noise-free version of the 
hybrid stochastic system (I13I14I ). 

Let us now examine the contraction properties of equations 
^ and (Ull. 

We have first 



dy dy 



VL^V^VLV^ 



(37" - 37 + 1)I4 



Second, 



so that A,,ax(^^^) =372-37+l< l(forO<7< 1). 



9g, 



v^v^ - V 

oy ox 



df_ 



(xi) 








€(X3) 



Now observe that Amax {^) g — I — —y^ < 1. Since V 
is an orthonormal projection, one then has Amax (^^) — 
Therefore, if 

37^ - 37 + 1 < e~2^ (15) 
then theorem |4] together with the corollaries of section IIV-BI 



imply that, after exponential transients, 

E(l|yir) < 



r.n2; / 27^^ + (1 - /?)(! + /3 - Pe^-)e^-a^, 



2(l-/3)(l-/3e2-) 
where /5 = 87^ - 87 + 1. 
To conclude, observe that 

3 



|Vx|p = i^||Rx,+i-x, 



Define the phase-locking quality 6 by 

3 

,5 = ^||Rx,+i-x,|i 



then one finally obtains 

E{S) < 



672^2 + 8(1 - + /? - /3e2^)e2v2 



2(l-/3)(l-/3e2-) 
after exponential transients. 

A numerical simulation is provided in Fig. 1. 



(16) 





Fig. 1. Numerical simulation using the Euler-Maruyama algorithm [15]. The 
following set of parameters was used: Cc = 0.1, ad = 0.05, t = 0.1. 
Two coupling strengths were tested: 7woak = 0.01 for plots (a), (b), (c), 
and 7strong = 0.2 for plots (c), (d), (e). Note that 7wcak does not satisfy 
condition )15) , while 7strong does, and yields the theoretical bound ~ 0.446 
(as provided by )16t ) on the phase-locking quality S. Plots (a) and (d) show 
the 2d trace of sample trajectories of the three oscillators for t G [0, 1]. Plots 
(b) and (e) show sample trajectories of the first coordinates of xi, R(x2) and 
(x3 ) as functions of time. Plot (c) and (f) show three sample trajectories 
of 5. 



VI. Perspectives 

We are now focusing on the following directions of re- 
search: 

• proving the state-dependent-metrics version of the con- 
tinuous and hybrid stochastic contraction theorems, 

• developping more elaborate conditions on dwell-times, 
and also hybrid switched versions of the theorems, 

• applying the synchronization-by-discrete-couplings anal- 
ysis to other types of coupled dynamical systems, 

• studying the robustness of hybrid controllers and ob- 
servers against random perturbations (for instance, the 
discrete observer for inertial navigation developped 
in [16]). 
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